	/*
	global directory1 "/Users/leandrocarvalho/Documents/UKB"
	global directory2 "/Users/leandrocarvalho/Dropbox/UK Biobank/paper"
	*/
	
	global directory1 "/Users/lcarvalh/Documents/UKB"
	global directory2 "/Users/lcarvalh/Dropbox/UK Biobank/paper"
	*/

	*global bw 			// 3652
	global polynomial 	"linear" // quadratic	
	global controls 	"male age age2 mixed_ethnicity asian black other_ethnicity NR_ethnicity Wales Scotland" 
	global MoB 	"Jan Feb Mar Apr May Jun Jul Aug Oct Nov Dec" 
	
	***************************************
	* 				Trends 				  *
	***************************************
	
	if "$polynomial" == "linear" {
		global posttrends	"c.DoBafter" 	
		global pretrends	"c.DoB"
	}
	
	if "$polynomial" == "quadratic" {
		global posttrends	"c.DoBafter##c.DoBafter" 	
		global pretrends	"c.DoB##c.DoB"
	}

	***************************
	* RD graphs
	***************************

	cap program drop RD_graph
	program define RD_graph
	
		preserve
			*1: dependent variable
			*2: title graph
			*3: sample restriction
			*4: name of file
			*5: y-axis title
			
			`3'
			
			qui {
	
				drop DoB
				gen temp = (year_birth*12 + month_birth) - (1957*12 + 9)
				gen DoB = floor(temp/3)
				
				reg `1' $pretrends	if DoB <  0 [aw=w]
					gen trend0 = _b[_cons] + _b[DoB]*DoB 								if DoB <= 0
				reg `1' $pretrends  if DoB >= 0 [aw=w]
					gen trend1 = _b[_cons] + _b[DoB]*DoB 								if DoB >= 0
	
								
				gen pop = 1
				collapse (mean) `1' (mean) trend0 (mean) trend1 (sum) pop, by(DoB)

			}	
				twoway 	(scatter `1' DoB [fw = pop], msymbol(Oh)	mcolor(black)) ///
						(line  trend0 DoB if DoB <= 0  , lcolor(black))  ///
						(line  trend1 DoB if DoB >= 0  , lcolor(black))  ///
						, ytitle("`5'") legend(off) ///
						xtitle("Quarter of Birth") xline(0,  lpattern(dash) lcolor(gs8)) ///
						graphregion(fcolor(white)) plotregion(lwidth(none)) scheme(s1mono) ylabel(, nogrid) /// 
						title("`2'", size(medium)) xscale(range(-20 19)) xlabel(-20(4)19) 
		
		restore	
	
	end		
	
	***************************
	* Tables
	***************************

	cap program drop tables
	program define tables
	
		*1: reduced-form (RF) or 2SLS
		*2: list of outcome variables
		*3: name of file
	
		foreach y in `2' {
			
			mean `y' [aw=w] if DoB >= -365 & DoB <= -1
			
			if "`1'" == "RF" {
			
				regress 		`y' after 	$pretrends $posttrends $MoB	 							[aw=w], robust

				regress 		`y' after 	$pretrends $posttrends $MoB $controls 					[aw=w], robust

			}
			if "`1'" == "2SLS" {
			
				ivregress 2sls	`y' 		$pretrends $posttrends $MoB	 		  (edu16 = after)	[aw=w], robust

				ivregress 2sls	`y' 		$pretrends $posttrends $MoB	$controls (edu16 = after)	[aw=w], robust

			}
		
		}	
	
	end
	
	***************************
	* CDFs Graphs
	***************************
	
	cap program drop CDFs
	program define CDFs
	

		preserve
			*1: sample restriction
			*2: dependent variable
			*3: grid
			*4: if vertical lines needed
			*5: title x-axis
			*6: graph title

			qui {
			
				`1'
				drop if `2' == .
				
				gen x=.
				gen F0=.
				gen F1=.
			
				// Estimating difference in cumulative distributions
				local j = 1
				foreach i of numlist `3' {	
				
					replace x=`i' in `j'
					gen tempMeasure=(`2' <= `i') if `2' ~=.

					regress tempMeasure $pretrends if edu16 == 0 & DoB < 0 [aw=w], robust // 
						replace F0=_b[_cons] 					in `j'
						local F0 = _b[_cons]
					ivregress 2sls tempMeasure $pretrends $posttrends $MoB (edu16 = after) [aw=w], robust //  
						replace F1=`F0' + _b[edu16]				in `j'
						
					drop tempMeasure
					
					local j = `j' + 1
				}
			}
			
			twoway  (line F0 x, lcolor(black) lpattern(shortdash)) || ///
					(line F1 x, lcolor(red)) || ///
					, `4' scheme(s1mono) graphregion(fcolor(white)) title("", size(medium))    ///
					xtitle("`5'", size(medium))  /// 
					ylabel(, nogrid) legend(order(1 "Pre" 2 "Post") row(2) ring(0) position(5) region(lwidth(none)))

		restore
			
	end

	***************************
	* Figure 1
	***************************	
	
	RD_graph edu16 				""		""	"Figure 1" "Fraction Staying in School until Age 16"
	
	***************************
	* Table 1
	***************************	
		
	tables RF 					"edu16 no_qualification CSE Olevel Alevel college" 	"Table 1"

	***************************
	* Figure 2
	***************************	
	
	RD_graph CSE_Olevel 		""	""	"Figure 2" "Fraction with a CSE or O-level Qualification"

	***************************
	* Table 2
	***************************	
	
	tables RF 					"index_anthro index_spiro index_bp index_health" 	"Table2 RF"	
	tables 2SLS 				"index_anthro index_spiro index_bp index_health" 	"Table2 2SLS"	
	
	***************************
	* Figure 3
	***************************
			
	RD_graph index_anthro  		"Body Size Index" 	""	"Figure 3A" ""
	RD_graph index_spiro  		"Lung Function Index" 			""	"Figure 3B" ""
	RD_graph index_bp  			"Blood Pressure Index" 		""	"Figure 3C" ""
	RD_graph index_health  		"Summary Index" 			""	"Figure 3D" ""
		
	***************************
	* Figure 4
	***************************

	gen BMI30 = (raw_BMI <= 30) if raw_BMI < .
	
	RD_graph BMI30 		""	""	"Figure 4" "Fraction with a BMI Below 30"
	
	***************************
	* Figure 5
	***************************

	_pctile raw_BMI [aw=w] if DoB < 0, p(1 99)
	local min = r(r1)
	local max = r(r2)
	local delta = (`max' - `min') / 20
	
	CDFs ""	raw_BMI		"`min'(`delta')`max'"	"xline(25, lwidth(vthin) lcolor(gs8)) xline(30, lwidth(vthin) lcolor(gs8))"	"Body Mass Index"	"Figure 7" 
	
	***************************
	* Figures 6-8
	***************************

	global min = 0
	global max = 0
	foreach y in index_anthro index_bp index_spiro { // index_health {
		_pctile `y' [aw=w] if DoB < 0, p(1 99)
		global min = min($min,r(r1))
		global max = max($max,r(r2))
	}
	global delta = ($max - $min) / 20
		
			*1: sample restriction
			*2: dependent variable
			*3: grid
			*4: if vertical lines needed
			*5: file name of graph #1
			*6: graph title
		
	CDFs ""	index_anthro	"$min($delta)$max"	""	"Body Size Index"	"Figure 4" 
	CDFs ""	index_spiro		"$min($delta)$max"	""	"Lung Function Index"		"Figure 5" 
	CDFs ""	index_bp		"$min($delta)$max"	""	"Blood Pressure Index"	"Figure 6" 
			
	***************************
	* Table 3
	***************************
	
	cap program drop hypothesis_test
	program define hypothesis_test 
	
		*1: outcome variable
		*2: outcome name
	
		dis "*** `2'"
		
		preserve
		
			qui drop if `1' >= .
					
			qui WAD_ROSLA_RF_gen_lin `1' DoB $MoB, g(test) bw(1825)
			qui sum test
			local thetest=r(mean)
			qui use "data/secondary/Wad_stats_5yr_lin_t1957_`2'_full.dta", clear
			qui gen full_distr=(thestats>`thetest')
			sum full_distr			
			

		restore		
		
		preserve
		
			qui drop if `1' >= .
					
			qui WAD_ROSLA_RF_gen_lin `1' DoB $MoB, g(test) bw(1825) pctu(50)
			qui sum test
			local thetest=r(mean)
			qui use "data/secondary/Wad_stats_5yr_lin_t1957_`2'_half.dta", clear
			qui gen lower_half=(thestats>`thetest')
			sum lower_half
			
		restore		
		
		preserve
		
			qui drop if `1' >= .
					
			qui WAD_ROSLA_RF_gen_lin `1' DoB $MoB, g(test) bw(1825) pctl(50)
			qui sum test
			local thetest=r(mean)
			qui use "data/secondary/Wad_stats_5yr_lin_t1957_`2'_half.dta", clear
			qui gen upper_half=(thestats>`thetest')
			sum upper_half
						
		restore		
		
	end

		hypothesis_test index_anthro body_size
		hypothesis_test index_spiro  lung_function
		hypothesis_test index_bp     blood_pressure
		
			
	***************************
	* Figure 9
	***************************

	_pctile raw_bpd [aw=w] if DoB < 0, p(1 99)
	local min = r(r1)
	local max = r(r2)
	local delta = (`max' - `min') / 20

	CDFs ""	raw_bpd		"`min'(`delta')`max'"	"xline(80, lwidth(vthin) lcolor(gs8)) xline(90, lwidth(vthin) lcolor(gs8))"	"Blood Pressure Diastolic"	"Figure 8" 
	
